Basic_LOC01_plots

Basic LOC01 plots from Jennie’s python notebook, translated to R.

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.4.4     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.0
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(skimr)
library(patchwork)
#library(maps)
library(janitor)

Attaching package: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test

library(lubridate)
library(leaflet)
library(mapdeck)

Attaching package: 'mapdeck'

The following object is masked from 'package:tibble':

    add_column

theme_set(theme_classic())

Data

loc01_underway <- "data/LOC-01_Underway_continuous.txt"
loc01_drifter_comp <- "data/drifters_interp.csv"
loc01_ctd <- "data/CTD_downcast_upcast.csv"

Underway

uw_df <- read_csv(loc01_underway,
                     na = c("", "NA", "NaN"),
                     col_types = "cnninnnnnncni",
                     ) |> 
  mutate(DateTime_UTC = as.POSIXct(DateTime_UTC, 
                                   format = "%d-%b-%Y %H:%M:%S", 
                                   tz = "UTC"),
         Salinity_PSU = ifelse(Salinity_PSU < 29, NA, Salinity_PSU))


uw_dye_df <- uw_df |> 
  filter(DateTime_UTC > as.POSIXct("2023-09-02 12:00:00", tz = "UTC"),
         DateTime_UTC < as.POSIXct("2023-09-04 01:30:00", tz = "UTC"))

CTD

ctd_df <- read_csv(loc01_ctd) |> 
  clean_names()

Rows: 144400 Columns: 20
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr   (4): InOut, Upcast/Downcast, Start Date, flag:  0.000e+00
dbl  (15): Station, Cast, bpos, longitude [deg], Longitude [deg], t090C [ITS...
time  (1): Start Time UTC

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Need per-cast info: station, cast, in-out, start timestamp, start lat, start lon

ctd_loc <- ctd_df |> 
  mutate(timestamp = as.POSIXct(paste(start_date, as.character(start_time_utc)), format="%m/%d/%y %H:%M:%S")) |> 
  select(station, cast, in_out, 
         timestamp,
         lat = longitude_deg, lon = longitude_deg_2) |> 
  distinct(station, cast, .keep_all = TRUE)
ctd_loc

# A tibble: 14 × 6
   station  cast in_out timestamp             lat   lon
     <dbl> <dbl> <chr>  <dttm>              <dbl> <dbl>
 1       1     1 Out    2023-09-02 09:23:55  41.1 -70.7
 2       1     2 Out    2023-09-02 10:15:23  41.1 -70.7
 3       1     3 Out    2023-09-02 11:07:07  41.1 -70.7
 4       2     2 In     2023-09-02 14:12:08  41.1 -70.7
 5       2     1 Out    2023-09-02 14:50:52  41.1 -70.7
 6       3     1 In     2023-09-02 18:22:37  41.1 -70.7
 7       4     1 In     2023-09-02 23:53:32  41.1 -70.7
 8       5     1 In     2023-09-03 04:34:53  41.1 -70.7
 9       6     1 In     2023-09-03 04:53:02  41.1 -70.7
10       7     1 In     2023-09-03 05:13:51  41.1 -70.7
11       8     1 In     2023-09-03 09:41:50  41.1 -70.7
12       9     2 In     2023-09-03 13:28:13  41.1 -70.7
13      10     1 In     2023-09-03 18:29:36  41.1 -70.7
14      11     1 In     2023-09-03 23:16:45  41.0 -70.7

#|fig-height: 800
ggplot(ctd_df, aes(dep_sm_salt_water_m, rhodfl_tc0_ppb, color = upcast_downcast)) +
  geom_line() +
  coord_flip() +
  scale_y_log10() +
  scale_x_reverse() +
  facet_grid(rows = vars(in_out), cols = vars(station, cast))

Drifters

drifter_df <- read_csv(loc01_drifter_comp)

Rows: 2299 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): asset
dbl  (6): lat, lon, deployment, dye, temp, sal
lgl  (1): interpolated
dttm (1): datetime

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

drifter_df |> 
  ggplot(aes(datetime, dye, color = dye)) +
  geom_point() +
  scale_y_log10()

Warning in self$trans$transform(x): NaNs produced

Warning: Transformation introduced infinite values in continuous y-axis

Warning: Removed 1 rows containing missing values (`geom_point()`).

# Create a palette that maps factor levels to colors
pal_fac <- colorFactor(c("red", "navy"), domain = ctd_loc$in_out)

#pal_rho <- colorQuantile("magma", drifter_df[["Dye_ppb"]], n = 10)
pal_rho <- colorQuantile("magma", uw_dye_df[["sensor_ppb"]], n = 50)

leaflet(drifter_df) |> 
  addTiles() |> 
  addCircleMarkers(lng = ~lon,
             lat = ~lat,
             color = ~pal_rho(dye),
             popup = ~asset,
             radius = 2,
             stroke = FALSE,
             fillOpacity = 0.5) |> 
  addCircleMarkers(
    data = uw_dye_df,
    radius = 2,
    stroke = FALSE,
    fillOpacity = 0.5,
    color = ~pal_rho(uw_dye_df[["Dye_ppb"]])) |> 
  addCircleMarkers(
    data = ctd_loc,
    color = ~pal_fac(in_out)
  )

Assuming "lon" and "lat" are longitude and latitude, respectively

Assuming "Longitude" and "Latitude" are longitude and latitude, respectively

Timeseries

Underway

dye_plot <- ggplot(uw_dye_df, aes(DateTime_UTC, Dye_ppb)) +
  geom_line() +
  scale_y_log10() +
  labs(x = "Time",
       y = "Dye (ppb)")

dye_plot_ctd <- ggplot(uw_dye_df, aes(DateTime_UTC, Dye_ppb)) +
  geom_point() +
  geom_line() +
  scale_y_log10() +
  labs(x = "Time",
       y = "Dye (ppb)")

temp_plot <- ggplot(uw_dye_df, aes(DateTime_UTC, Temp1_C)) +
  geom_line() +
  labs(x = "Time",
       y = "Temp (C)")
sal_plot <- ggplot(uw_dye_df, aes(DateTime_UTC, Salinity_PSU)) +
  geom_line() +
  labs(x = "Time",
       y = "Salinity (PSU)")
ta_plot <- ggplot(uw_dye_df, aes(DateTime_UTC, Corrected_TA_umol_kg_)) +
  geom_point(shape = 1) +
  labs(x = "Time",
       y = expression(paste(TA ~ (μmol ~ kg^-1))))

(dye_plot + temp_plot) / (sal_plot + ta_plot)

Warning: Removed 67280 rows containing missing values (`geom_point()`).

Maps

Underway

# Outline map
map_data = map_data("usa")

head(uw_dye_df)

# A tibble: 6 × 17
  DateTime_UTC        Latitude Longitude location_data_flag Dye_ppb    mv
  <dttm>                 <dbl>     <dbl>              <int>   <dbl> <dbl>
1 2023-09-02 12:00:02     41.1     -70.7                  1    1.75   872
2 2023-09-02 12:00:04     41.1     -70.7                  1    1.77   883
3 2023-09-02 12:00:06     41.1     -70.7                  1    1.71   852
4 2023-09-02 12:00:08     41.1     -70.7                  1    1.71   855
5 2023-09-02 12:00:10     41.1     -70.7                  1    1.73   862
6 2023-09-02 12:00:12     41.1     -70.7                  1    1.72   858
# ℹ 11 more variables: Temp1_C <dbl>, c1 <dbl>, Temp2_C <dbl>,
#   Salinity_PSU <dbl>, Location_File <chr>, Corrected_TA_umol_kg_ <dbl>,
#   Data_Flag <int>, fCO2_SW_SST_uatm <dbl>, fCO2_ATM_interpolated_uatm <dbl>,
#   dfCO2_uatm <dbl>, WOCE_QC_FLAG <dbl>

myscaler <- function(x, from, to) {
  high=1.2
  low=0.1
  ifelse(x<low,0,ifelse(x>high,1,(x-low)/(high-low)))
}


# Plot 
ggplot() + 
  geom_map(data = map_data, map = map_data,
           aes(long, lat, map_id = region),
           fill = "white", color = "black") +
  geom_point(data = uw_dye_df, 
             aes(x = Longitude, y = Latitude, color = log10(Dye_ppb)),
             size = .4) +
  scale_color_viridis_c(option = "D",
                        name = "Dye (log ppb)",
                        #trans = scales::log10_trans(),
                        rescaler = myscaler) +
  xlim(-70.8, -70.6) +
  ylim(41.01, 41.15) +
  coord_quickmap() +
  ggtitle("Dye") +
  xlab("Longitude") +
  ylab("Latitude")

Warning in geom_map(data = map_data, map = map_data, aes(long, lat, map_id =
region), : Ignoring unknown aesthetics: x and y

ggplot() + 
  geom_map(data = map_data, map = map_data,
           aes(long, lat, map_id = region),
           fill = "white", color = "black") +
  geom_point(data = uw_dye_df, 
             aes(x = Longitude, y = Latitude, color = Salinity_PSU),
             size = .4) +
  scale_color_viridis_c(option = "D") +
  xlim(-70.8, -70.6) +
  ylim(41.01, 41.15) +
  coord_quickmap() +
  ggtitle("Salinity") +
  xlab("Longitude") +
  ylab("Latitude")

Warning in geom_map(data = map_data, map = map_data, aes(long, lat, map_id =
region), : Ignoring unknown aesthetics: x and y

ggplot() + 
  geom_map(data = map_data, map = map_data,
           aes(long, lat, map_id = region),
           fill = "white", color = "black") +
  geom_point(data = uw_dye_df, 
             aes(x = Longitude, y = Latitude, color = Temp1_C),
             size = .4) +
  scale_color_viridis_c(option = "B",
                        name = "Temp (C)") +
  xlim(-70.8, -70.6) +
  ylim(41.01, 41.15) +
  coord_quickmap() +
  ggtitle("Temperature") +
  xlab("Longitude") +
  ylab("Latitude")

Warning in geom_map(data = map_data, map = map_data, aes(long, lat, map_id =
region), : Ignoring unknown aesthetics: x and y

cruise_TA <- uw_df |> 
  filter(!is.na(Corrected_TA_umol_kg_)) |> 
ggplot() + 
  geom_map(data = map_data, map = map_data,
           aes(long, lat, map_id = region),
           fill = "white", color = "black") +
  geom_point(aes(x = Longitude, y = Latitude, color = Corrected_TA_umol_kg_),
             size = 2) +
  scale_color_viridis_c(option = "D",
                        name = expression(paste(TA ~ (μmol ~ kg^-1)))) +
  xlim(-72.1, -70.4) +
  ylim(40.5, 41.5) +
  coord_quickmap() +
  ggtitle("Alkalinity") +
  xlab("Longitude") +
  ylab("Latitude")

Warning in geom_map(data = map_data, map = map_data, aes(long, lat, map_id =
region), : Ignoring unknown aesthetics: x and y

dye_TA <- uw_dye_df |>  
  filter(!is.na(Corrected_TA_umol_kg_)) |> 
  ggplot() +
  geom_map(data = map_data, map = map_data,
           aes(long, lat, map_id = region),
           fill = "white", color = "black") +
  geom_point(aes(x = Longitude, y = Latitude, color = Corrected_TA_umol_kg_),
             size = 2, shape = 1) +
  scale_color_viridis_c(option = "D",
                        name = expression(paste(TA ~ (μmol ~ kg^-1)))) +
  xlim(-70.8, -70.6) +
  ylim(41.01, 41.15) +
  coord_quickmap() +
  ggtitle("Alkalinity") +
  xlab("Longitude") +
  ylab("Latitude")

Warning in geom_map(data = map_data, map = map_data, aes(long, lat, map_id =
region), : Ignoring unknown aesthetics: x and y

cruise_TA + dye_TA

Drifters

  ggplot() +
  geom_point(data = uw_dye_df, 
             aes(x = Longitude, y = Latitude, color = log10(Dye_ppb)),
             size = .4) +
  geom_point(data = drifter_df, aes(lon, lat, shape = asset, color = log10(dye))) +
  scale_color_viridis_c(option = "D",
                        name = "Dye (log ppb)",
                        trans = scales::log10_trans()) +
  geom_point() +
  coord_quickmap()

Warning in FUN(X[[i]], ...): NaNs produced

Warning in self$trans$transform(x): NaNs produced

Warning: Transformation introduced infinite values in discrete y-axis

drifter_df |> 
  ggplot(aes(datetime, dye)) +
  geom_line(data = uw_dye_df, aes(DateTime_UTC, Dye_ppb), color = "lightblue") +
  geom_line() +
  scale_y_log10() +
  facet_grid(rows = vars(asset))

Warning in self$trans$transform(x): NaNs produced

Warning: Transformation introduced infinite values in continuous y-axis

Mapdeck map

#m <- grDevices::colorRamp(c("blue", "white", "yellow"), bias = 6)(df[["Dye_ppb"]])
#pal <- col_quantile("viridis", df[["Dye_ppb"]], n = 30)(df[["Dye_ppb"]])

uw <- uw_dye_df |> 
  mutate(dye = log10(Dye_ppb))
         
dr <- drifter_df |> 
  mutate(dye = log10(dye))

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `dye = log10(dye)`.
Caused by warning:
! NaNs produced

key <- 'pk.eyJ1IjoiYmxvbmd3b3J0aCIsImEiOiJjbHI3NmtzMmYyZTBuMmluOGZ4dmhuaWFqIn0.eW7yUb_wqL7LQJOURH9IRQ'
mapdeck(token = key, style = mapdeck_style("light")) |> 
  add_scatterplot(
    data = uw, 
    lat = "Latitude", 
    lon = "Longitude",
    radius = 10, 
    fill_colour = "dye",
    layer_id = "dye_layer",
    palette = "viridis") |> 
  add_scatterplot(
    data = dr, 
    lat = "lat", 
    lon = "lon",
    radius = 10, 
    fill_colour = "dye",
    layer_id = "drifter_layer",
    palette = "viridis")

Registered S3 method overwritten by 'jsonify':
  method     from    
  print.json jsonlite